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Abstract 

We analyze the density distribution and the adsorption of solvent hard spheres in an annular 
slit formed by two large solute spheres or a large solute and a wall at close distances by means 
of fundamental measure density functional theory, anisotropic integral equations and simulations. 
We find that the main features of the density distribution in the slit are described by an effective, 
two-dimensional system of disks in the vicinity of a central obstacle. This has an immediate 
consequence for the depletion force between the solutes (or the wall and the solute), since the 
latter receives a strong line-tension contribution due to the adsorption of the effective disks at the 
circumference of the central obstacle. For large solute-solvent size ratios, the resulting depletion 
force has a straightforward geometrical interpretation which gives a precise "colloidal" limit for 
the depletion interaction. For intermediate size ratios 5 ... 10 and high solvent packing fractions 
larger than 0.4, the explicit density functional results show a deep attractive well for the depletion 
potential at solute contact, possibly indicating demixing in a binary mixture at low solute and high 
solvent packing fraction. 
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I. INTRODUCTION 



The equilibrium statistical theory of inhomogeneous fluids whose foundations were laid 
out in the 1960's is a well-studied subject which has been developed since by a fruitful 
interplay between simulations and theory. For the simplest molecular model, the hard-body 
fluid mixture, this has led to the development of powerful, geometry-based density function- 
als 5 la, lZI, |8| which accurately describe adsorption phenomena, phase transitions (such 

as freezing for pure hard spheres or demixing for entropic colloid-polymer mixtures) and 
molecular layering near obstacles. Such a class of density functionals has not been found yet 
for fluids with attractions, however, significant progess has been achieved in the description 
of bulk correlation functions and phase diagrams through the method of integral equations 



Strong inhomogeneities occur if fluids are confined on molecular scales, a topic which 
enjoys continuous interest jl2|, [l^. For example, the packing of molecules between parallel 
walls leads to oscillating forces between them which can be measured experimentally and 
determined theoretically [l^ . With the development of preparational techniques for colloids 
and their mixtures, the "molecular" length scale has been lifted to the range of several nm 
up to fim which even allows the direct observation of modulated density profiles through 
microscopy besides the measurement of resulting forces on the walls. Furthermore, in the 
colloidal domain the possibility to tailor the interparticle interactions to a certain degree 
allows the close realization of some of the favourite models for simple fluids fancied by theo- 

n 

rists, such as e.g. hard spheres [15|. This has opened the route to quantitative comparisons 
between experiment, theory and simulations. 

In an asymmetric colloidal mixture with small "solvent" particles and at least one species 
of larger solute particles, the phenomenon of solvent-mediated, effective interactions between 
the solute particles may give rise to various phase separation phenomena 
theoretical point of view, these effective interactions are interesting since they facilitate 
the description of mixtures in terms of an equivalent theory for one species interacting by 
an effective potential fl?! . If the solvent particles possess hard (or at least steeply 

repulsive) cores, they are excluded from the region between two solute particles if the latter 
are separated by less than one solvent diameter. This gives rise to strong depletion forces 
whose understanding is crucial for the concept of an effective theory containing only solute 



16|. From a 
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FIG. 1: (color online) The annular wedge which is formed between two large solute particles for 
separations h < a. Black areas denote domains which are forbidden for the centers of the solvent 
spheres. 

degrees of freedom. Since the magnitude of the depletion force is directly linked to the 
solvent density distribution around the solutes, the quantitative investigation of the solvent 
confined between solute particles appears to be important. The confinement becomes rather 
extreme for large asymmetry between solute and solvent (see below). 

In the following we want to concentrate on the idealized system of additive hard spheres 
and in particular on the effective interaction between two solute particles in solvent, i.e. the 
case of infinite dilution of solute particles in the colloidal mixture. The case of the interac- 
tion of one solute particle with a wall is a special case in that the radius of the other solute 
particle goes to infinity. For the case of hard spheres (solvent diameter a = 2Rs, solute ra- 
dius Rb such that the size ratio is a = Rb/Rs and R = Rb + Rs is the radius of the exclusion 
sphere for solvent centres around a solute sphere) it has been found that the theoretical 
description using "bulk" methods becomes increasingly inaccurate for size ratios a > 5 and 
solvent densities p* = psU^ ^ 0.6, when compared to simulations. Here, the terminus "bulk" 
methods refers to methods which determine the bulk pair correlation function between the 
colloids, 5'bfe(r), from which the depletion potential is obtained as (3W = — In Qbb- Bulk inte- 
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IE) methods [19|] and the "insertion" trick in density functional theory (DFT) 



231] fall into this category. For larger size ratio, one would expect that the well- 



known Derjaguin approximation [2J, |25|, |26| becomes accurate very quickly. Interestingly, 



however, both the mentioned bulk methods and simulations deviate significantly from the 
Derjaguin approximation (besides disagreeing with each other) for size ratios 10 and solvent 
densities p* = 0.6 ... 0.7 26|, |27|, |28|| when the surface-to-surface separation h between the 
solute particles is close to a {h < a, i.e. at the onset of the depletion region where the solvent 
particles are "squeezed out" between the solute particles). One notices that ioi h < a the 
solute particles form an annular wedge with a sharp edge which restricts the solvent particles 
to quasi-2d motion (see Fig. [T]). Taking this observation into account, the phenomenological 



analysis of Ref. [26j predicted that the solvent adsorption at this edge leads to a line contri- 
bution to the depletion potential which is proportional to the circumference of the circular 
egde and thus to i?^/^. Such a term in the depletion potential causes a very slow approach 
to the Derjaguin limit for large solutes (the Derjaguin approximation describes the deple- 



tion potential essentially 
pertaining to the solutes 



)y volume and area terms of the overlap of the exclusion spheres 



25 



26| . see below). In this way, a generalized Derjaguin approx- 



imation serving as a new "colloidal" limit can be formulated which turns out to have a far 
more general meaning 29|| than anticipated. Using the concept of morphological (morpho- 
metric) (morphometric) (morphometric) (morphometric) (morphometric) (morphometric) 
(morphometric) (morphometric) thermodynamics introduced in Ref. 30|], one finds that the 
insertion free energy of two solute particles (and thus the depletion potential between them) 
only depends on the volume, surface area, and the integrated mean and Gaussian curvatures 
of the solvent accessible surface around the two solutes. The coefficients of these four terms 
depend on the solvent density hut not on the specific type of surface. In this manner, the 
phenomenological line tension of Ref. 2^ is related to the general coefficient of mean cur- 
vature of the hard-sphere fluid [29|, and the Derjaguin approximation is equivalent to the 
morphometric analysis restricted to volume and surface area terms. 

Since the depletion force in the hard-sphere system is directly linked to an integral over 
the solvent contact density on one solute (see Eqs. ([1]) and ([2]) below), one should be able to 
connect the morphometric approach to features of the solvent density proflle in the annular 
wedge. This is a strong motivation for us to explicitly determine the wedge density proflle. 
We will do so by means of density functional theory and anisotropic integral equations, and 
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compare it to results of Monte-Carlo (MC) simulations for selected parameters. For the 
intermediate densities p* = 0.6 and 0.7 and various size ratios between solute and solvent 
results for the depletion force from such explicit DFT calculations were already presented 



in Ref. 



29|], and good agreement with the morphometric depletion force was obtained for 



size ratios between 5 and 40. Here, we will present a detailed analysis of the wedge density 
profiles and consider also higher densities. We will give a thorough comparison to recent 
results from MC simulations which have been obtained for the wall-solute interaction at 
a solvent density p* = 0.764 (r/s = {Ti/Q)psa^ = 0.4) and size ratios between 10 and 100 
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The paper is structured as follows. In Sec. [TTl we introduce the basic notions of density 
functional and integral equation theory and present a short description of the Monte Carlo 
simulations employed here. In Sec. IIII Al we analyze in detail the density profiles and the 
corresponding depletion forces for the wall-sphere geometry for the particular solvent density 
p* = 0.764. This geometry permits explicit calculations up to solute-solvent size ratios 100 
and a test of the "colloidal limit" of the morphometric depletion force. Sec. IIIIBI gives an 
analysis of the depletion force and potential in the sphere-sphere geometry (corresponding 
to the important case of the effective interaction in a dilute mixture) for intermediate size 
ratios 5 and 10 and higher solvent densities p* = 0.8 . . . 0.9. 



II. THEORY AND METHODS 



We consider two scenarios (see Fig. [2]): (a) One hard solute sphere immersed in the hard 
sphere solvent of density ps confined to a slit, created by two hard walls at distance L. 
The surface-to-surface distance between one wall and the solute sphere is denoted by h, 
furthermore L ^ h such that the correlations from one wall do not influence the solute 
interaction with the other wall, (b) Two hard solute spheres immersed in the bulk solvent 
spheres with surface-to-surface distance h. The solvent accessible surface is given by the 
dashed lines in Fig. [21 thereby one sees that for h < a the two solutes (or the solute and 
one wall) form an annular wedge in which the solvent adsorbs. The solvent density profile 
p(r) = p(r||, z) depends only on z, the coordinate on the symmetry axis and the distance r\\ 
to the symmetry axis. In the latter case (b), when the first colloid is centered at the origin 
and the second one aX z = 2/?^ + /i, the total force f{h) on the first colloid is obtained by 
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integrating the force —Vubsir) {uhs is the sohite-solvent potential) between the solute and 
one solvent sphere over the density distribution p(r): 

(3f{h) = - Jd'r p(r) e, ■ Vn,,(r) (1) 

= 2ttR^ j d(cos^)cos^ p(r) , 
[|r| = R , cos 6' = r ■ e^] . 

Thus the total depletion force reduces to an integral over essentially the contact density on 
the exclusion sphere around the colloid. This follows from exp(— /5ubs)/5Vubs = — f(5(|r| —R) 
and the observation that p(r) exp(/5ub<() is continuous across the exclusion sphere surface. In 
case (a), the force on the colloid is the negative of the excess force on the wall and the 
excess force is determined by the total force on the wall minus the contribution from the 
bulk solvent pressure p. Through the wall theorem, the latter is given by (3p = where p^ 
is the contact density of the solvent at a single wall. By an argument similar to the above 
one and putting the exclusion surface of the wall at 2; = 0, is determined as 

/•oo 

/3fM = -2ti / r|| dr\\ {p{r\\,z = 0) - p^) . (2) 
Thus, fu) is equivalent to the excess adsorption at the wall which makes it somewhat easier 



to determine in MC simulations than / 



32|. 



A. Density functional theory 

The equilibrium solvent density profile p(r) = Pcq(i") can be determined directly from the 
basic equations of density functional theory. The grand potential functional is given by 

n[p] = P'^lp] + ^^1p] - J dr{p - V-*(r)) , (3) 

where JF^'^ and denote the ideal and excess free energy functionals of the solvent. The 
solvent chemical potential is denoted by p and the solute(s) and/or the walls define the 
external potential The ideal part of the free energy is given by 

/3r^ = J rfrp(r)(ln(p(r)A3)-l) , (4) 

with A denoting the de-Broglie wavelength. The equilibrium density profile Peq(r) for the 
solvent at chemical potential p = /?~^ ln(p5 A'^) + p"^^ (corresponding to the bulk density ps) 
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FIG. 2: View of the geometric configurations used in this work, (a) Solute sphere of radius Rb 
immersed in a solvent-filled slit of width L. Only the density profile in the annular wedge between 
the left wall and the solute is of interest since it determines the depletion force between solute and 
one wall. Note that one can also determine the slit density profile and the corresponding depletion 
force for 2i?b > L (i.e. when the solute does not fit into the slit) as long as L is large enough that 
the correlations from the right wall do not reach into the annular wedge, (b) Two solute spheres of 
radius Rb at distance h immersed in bulk solvent. For both setups, the solvent sphere diameter is 
given by a, and the radius of the exclusion sphere around a solute particle is given hy R = Rb + cr/2. 



is found by minimizing the grand potential in Eq. ([3]): 



ps 



5piv) 



/3p' 



(5) 



For an explicit solution, it is necessary to specify the excess part of the free energy. Here 
we employ two functionals of fundamental measure type (FMT). These are the original 
Rosenfeld functional (FMT-RF) and the White Bear functional (FMT-WB, with Mark 
I from Refs. 0, and Mark II from Ref. ^), for another closely related variant see Ref. [3]. 
It has been demonstrated that FMT gives very precise density profiles also for high densities 
of the hard sphere fluid in various circumstances. Furthermore, the White Bear II functional 



possesses a high degree of self-consistency with regard to scaled-particle considerations . 
However, we do not consider the tensor-weight modifications of these functionals which 
are necessary to obtain a correct description of the liquid-solid transition 3J] and are of 
higher consistency in confining situations which reduce the dimensionality of the system 
("dimensional crossover", see Ref. 
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36l|). This might be an issue in some circumstances 



(see below). 

Th explicit forms of the excess free energy are given in App. [Bl Numerically, due to the 
non-local nature of T'^^, Eq. ([5]) corresponds to an integral equation for the density profile 
depending on the two variables r\\ and z. The sharpness of the annular wedge in the solute- 
solute and the solute-wall problems for high size ratios necessitates rather fine gridding 
(which makes the calculation of J^°^ and 5J^^^/ 6p very time-consuming) and introduces an 
unusual slowing-down of standard iteration procedures of Picard type. To overcome these 
difficulties, we made use of fast Hankel transform techniques and more efficient iteration 
procedures. These are described in detail in App. [B] as well. 

We remark that there have been earlier attempts to obtain explicit density distributions 



using DFT around two fixed solutes and thereby to extract depletion forces 37|, 



38 



Ref. 



this was done using the simple Tarazona I functional for hard spheres [40|, and 



m. In 



considerations were limited to a small solvent packing fraction of 0.1 and solute-solvent size 



ratios of 5. In Ref. 



391], a minimization of FMT-RF was carried out using a real-space 



technique for solvent densities up to p* = 0.6 and size ratios 5. In that respect, the present 
technique is superior in that it allows us to present solutions for solvent densities up to 
p* = 0.9 and size ratios up to 100. 



B. Integral equations 

The central objects within the theory of integral equations are the correlation functions 
on the one and two-particle level in the solvent. We are interested in the two-particle 
correlation functions in the presence of an external background potential (one fixed object, 
solute or wall), from which the explicit solvent density profiles around the two fixed objects 
(solute-solute or solute-wall) follow (see below). This is usually referred to as the method of 
inhomogeneous or anisotropic integral equations, first employed for Lennard- Jones fiuids on 
solid substrates 4l|, |42] and for hard-spheres (HS) fiuids in contact with a single hard wall 
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in Ref. 43|]. Liquids confined between two parallel walls (a planar slit) have been extensively 



studied by Kjellander and Sarman 
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For an arbitrary background potential V{r) which corresponds to an equilibrium back- 
ground density profile for the solvent py(r), the pair correlation function 5fjj(r,ro) = 
/ijj(r, To) + 1 describes the normalized probability to find a particle of species i at posi- 
tion r if another particle of species j is fixed at position tq. For our system, the species 
index is either s (solvent particle) or b (big solute particle). In the solute-solute case, V is 
given by the potential of one solute particle, whereas in the solute-wall case it is the poten- 
tial of the wall. Then the equilibrium density profile discussed in the previous subsection 
is related to the pair correlation function through p(r) = py(r) (^^^(r, Tq) (fq specifies the 
position of the (other) solute particle). The depletion force follows then through Eqs. ([T]) 
and ^. 

The corresponding direct correlation functions of second order Qj(r, fq) are related to hij 
through the inhomogeneous Ornstein-Zernike (OZ) equations 



hij{r,ro) - Cij{r,Vo) = ^ J dr'pv,fc(r')^ife(r, r')cfcj(r', Fq) 

k=b,s 



(6) 



In the dilute limit for the solute particles which we consider here, the background density 
Pv^b for the solutes is zero, therefore the OZ equations reduce to 



hbs{r,ro) - Cbs{r,ro) 



dr'pvir')hssir, r')css{r', Tq) 



(7) 
(8) 



The background density profile is linked to Css through the Lovett-Mou-Buff-Wertheim 



equation 
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Vpy(r) = -ppv{r)\/Vir)+ pv{r) J rfr'c,,(r, r') Vpy(r') . 



(9) 



A third set of equations is necessary to close the system of equations. The diagrammatic 
analysis of Ref. |l| provides the general form of this closure which reads: 



liagij{r, Fq) + pUij{r - Tq) = hij{r, Fq) - Cij{r, Tq) - bij{r, Fq) 



(10) 



where the pair potentials in the solute-solvent mixture and bij denote the bridge 

functions specified by a certain class of diagrams which, however, can not be resummed in a 



closed form. For practical applications, these bridge functions need to be specified in terms 
of hij and Cij to arrive at a closed system of equations. Among the variety of empirical forms 
devised for this connection we mention those which, in our opinion, rely on somewhat more 
general arguments. These are the venerable hypernetted chain (HNC), Percus-Yevick (PY) 
closure and the mean spherical approximation (MSA) which can be derived by systematic 
diagrammatic arguments 48|, and the reference HNC closure which introduces a suitable 
reference system for obtaining bij with subsequent free energy minimization 49|. (For bulk 
properties of liquids, the self-consistent Ornstein-Zernike approximation (SCOZA) and 



the hierarchical reference theory (HRT) 10|] are very successful. They rely on a closure of 
MSA type but it appears to be difficult in generalizing them to inhomogeneous situations 
such as considered here.) In calculations, we considered the closures 



7i,-ln(l + 7.,) (PY), 

exp{[l - exp(-^ijr)]7ij} - 1 



lij - In 



1 + 



7^ 



1 - exp(-^ijr) 
(MV) . 



(RY) 



(11) 
(12) 

(13) 



2 1 + aij%j 

where = hij — Qj. The Percus-Yevick (PY) approximation Eq. ( ITT]) is exactly solvable in 



the bulk case even for the muticomponent HS fluid [50|. However, beyond providing good 
qualitative behavior PY does not produce precise quantitative results in general, since it fails 
at contact, where the value of pair distribution function is too small. For bulk systems, PY 
generates a noticeable thermodynamic inconsistency, i.e. a discrepancy between different 
routes to the equation of state. To overcome these deficiencies a refined approximation to 
the closure, which interpolates between HNC {bij = 0) and PY closures was suggested by 



Rogers and Young (RY) [Sli, see Eq. (fT2|) where the ^ij are adjustable parameters, which 



can be found from the thermodynamic consistency requirement 



5l| . An even more suitable 



for the asymmetric HS mixtures variant of the modified Verlet (MV) closure was suggested 



in Ref. 



52], see Eq. 0131) . Here the parameters tty are chosen to satisfy the exact relation 



between bij{0) and the third virial coefficient at low densities. For the choice of the auxiliary 
parameters C,ij and aij in the present work, see App. [C] 

We have solved the set of equations fl71)- f|T3l) for the solute-wall case where we could 
employ similar numerical methods as in the numerical treatment of density functional theory. 
(This is not possible for the solute-solute case, for an efficient method applicable in this case. 
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see Ref. |53|.) In the solute-wall case, the background density profile pi/(r) = Pv{z) depends 
only on the distance z from the wall. The two-particle correlation functions depend on the 
^-coordinates of the two particles individually and the difference in the radial coordinates: 
(7jj(r, To) = gij{z,Zo,r\^ — r|| o). Thus, the Ornstein-Zernike equations ([7j) and become 
matrix equations for the correlation functions in the ^-coordinates, and are diagonal for the 
Hankel transforms of the correlation functions in the r||-coordinates. For more details on 
the numerical procedure, see App. O 

C. Simulation Details 

The Monte Carlo simulations for the wall-solute system were performed at fixed particle 
number and volume of the simulation box. The upper boundary of the box was given by a 
large hard sphere, the lower boundary by a planar hard wall (see hatched regions in Fig. [3D. 
All remaining boundaries were treated as periodic. We sampled the configuration space of 
the small spheres by standard single particle translational moves. In order to impose the 
asymptotic solvent bulk density p* = 0.764 {tJs = 0.4), the concentration of small spheres 
in the box was set such that the density at contact with the hard wall far away from the 
wedge (i. e. in the region marked by grey squares in Fig. [3l averaged over a depth 0.02 
a) settled to pwaii = 4.88 within 1% error. This value pwaii = 4.88 was obtained from the 
density profile of hard spheres at a hard wall at the bulk packing fraction 0.4, calculated with 
FMT-WBII which is very accurate. System sizes ranged from 1800 particles for Rf, = 10 a 
to 8000 particles for Ri, = 25 a. In order to access the configurations inside the narrow part 
of the wedge with sufficient accuracy, large numbers of Monte Carlo sweeps were required. 
We equilibrated the systems for 5 x 10^ MC sweeps (i. e. attempted moves per particle) 
each. For data acquisition, we performed between 10^ and 10^ sweeps (depending on the 
parameters h and Rb) and averaged over 10^ to 10^ samples. Note that from the simulations 
we obtained only the wedge density profiles (see Figs. [5H7| below) but did not attempt to 
obtain the depletion force from Eq. ([2]) where one needs an excess adsorption integral at the 



wall. This would require considerably more sweeps 
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33| . see a. 



statistical errors of the simulated depletion force according to Ref. 



so Fig. [8] below for the 



32|. 
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FIG. 3: Sketch of the simulation box. Hard walls are marked by hatches. The remaining boundaries 
were treated as periodic. The grey squares mark the regions in which the density at wall contact 
was sampled. 



III. RESULTS 



A. Wall solute interaction 



In this section, we analyze case (a), the problem of a solute sphere close to a wall. Previous 



simulation work 
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331 ] has concentrated on the particular state point p* = 0.764 



(packing fraction 7]^ = 0.4), therefore we present explicit results for the density profiles also 
for this state point to facilitate comparison. 



1. Density profiles 

In Fig. m we show DFT results for the solvent density profile between wall and solute for 
a solute-solvent size ratio of 20 and solvent-wall distance h = (i.e. contact between wall 
and solute molecular surface; left panel) and h = a (i.e. the end of the depletion region; 
right panel). There is strong adsorption at the apex of the annular wedge, given by the 
coordinates z = and r|| = tq = \j2R{h — a) — {h — cr)^, as refiected by the main peak. 
Along the wall (z = 0), we observe strong structuring which is mainly dictated by packing 
considerations. The apex peak corresponds to a ring of solvent spheres centered at r|| = tq 
(for /i = 0) or just one sphere near r\\ = (for h = a). The second-highest peak in both 
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10. 




FIG. 4: (color online) Density profiles from FMT-WBII of solvent spheres confined to the annular 
wedge between a solute sphere (size ratio Rh/Rs = 20) and a wall. The bulk solvent density is 
p* = 0.764, and the wall separation L = 26 a (see Fig. [2] (a)). Grid resolution was equidistant in z 
with Az = 0.002 a and equidistant in x = ln(r||/cj) with Ax = 0.005 (see App. IB 2p . Left panel: 
solute-wall distance h = 0, right panel: h = a. 

panels of Fig. H] appears where the distance between wall and solute sphere, measured along 
the 2;-axis, is approximately 2o" and thus two rings of spheres fit between solute and wall. 
However, for /z = (left panel of Fig. H]) packing in r||-direction leads to further structuring 
of the density profile close to the apex of the wedge. 

The integral of the solvent density at the wall (i.e. the adsorption at the wall) directly 
determines the depletion force, see Eq. ([2]). Therefore, we analyze the density at the wall 
further. In Fig. [5] we show p{r\\,z = 0) computed with FMT-RF, FMT-WBII, lE-MV and 
compare the results to the simulation data of Ref. Overall, DFT performs quite well but 
it tends to overestimate the structuring of the profile. The currently most consistent version, 
FMT-WBII, improves over FMT-RF particularly in this respect. The integral equation 
results for the different closure are very similar to each other and are approximately of the 
same quality as the solutions of FMT-RF. Note that the simulation data of Ref. 3^ have 
been averaged over a distance of 0.02 a from the wall. This average was also applied to the 
DFT data (which are computed with mesh size Az = 0.002 a), whereas for the IE results we 
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FIG. 5: (color online) Comparison between the various theoretical methods and MC simulations 



(ours and from Ref. j3^) for the density at the wall (averaged over a distance of 0.02 a). Left 
panel: solute-wall distance h = 0, right panel: h = a. The solute-solvent ratio is 20 and the 
solvent bulk density is p* = 0.764. Note the logarithmic scale on the density axis. 

show a suitable weighted average of the zeroth and first bin {Az = 0.05), assuming linear 
interpolation. As seen in the right panel of Fig. [5l lE-MV gives a much lower density close 
to the apex of the wedge compared to the other methods. This is presumably due to the 
lower resolution in ^-direction which is possible for IE (see App. ICl) . 

To further investigate the issue of the dominant packing mode in the annular wedge, we 
performed DFT calculations also for the additonal size ratios 10, 30, 50 and 100 for the 
same solvent density p* = 0.764. Packing in r||-direction can be monitored conveniently 
by viewing the annular wedge close to the apex as a quasi-2d system (see Fig. [1] and the 
remarks in the Introduction). We define an effective, two-dimensional solvent density by 
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P2d(r||) 



dzp{rii,z) 



(14) 



where Zs{r^^) = R — {a — h) — \J — rj^ ^ rj^/(2_R) — [a — h) defines the surface of the 
exclusion sphere around the solute (and thus the radius-dependent width of the annular 
wedge). In Fig. [6] we show results for P2d('"||) for two configurations with slightly different 
solute-wall distances: h = 0.95 a (left panel) and h = a (right panel). In the latter case, 
one solvent sphere fits exactly between solute and wall at the apex, nevertheless p2d vanishes 
there due to the vanishing slit width, Zs{r\\ ^ 0) ^ 0. Note that this is also a consequence 
of an application of the potential distribution theorem 5^] which states that the 3d density 
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reaches a finite, maximum value of exp(/3yU^^) in a planar slit of vanishing width which is 
equivalent to a vanishing 2d density (//^^ is the excess chemical potential of the bulk solvent 
coupled to the slit). However, when the slit widens, the 2d density quickly rises 26|]. For 
h = a (right panel in Fig. [H]), this leads to the picture of an effective 2d system with a 
small, soft and repulsive obstacle centered at r|| =0 which induces moderate layering in the 
2d- "bulk" with an effective "bulk" density p2d,s ~ 0.7. . .0.8. For h = 0.95 cr (left panel in 
Fig. [6]), the obstacle in the center has a radial extent of tq ~ y^RaJTO which ranges from 
0.75 cr (size ratio 10) to 2.25 cr (size ratio 100) and induces stronger layering in the effective 
2d system. The oscillations occur around the same "bulk" density as in the case h = a. 
The simulation results show weaker oscillations than the DFT results. Since the first peak, 
e. g., occurs at wedge widths < 0.03 a, we are testing the dimensional crossover properties of 
DFT-FMT from 3d to 2d with a sensitive probe. The difference between the results of the 
simulations and the DFT versions employed here are in line with previous investigations on 



the dimensional crossover properties of FMT 



There it has been found that the strict 



2d limit of FMT-RF results in a somewhat peculiar (integrable) divergence in the hard disk 
direct correlation function C2d('") for r — and overestimated peaks in the corresponding 
structure factor. The tensor weight modifications introduced in Ref. result in better 
dimensional crossover properties j56| and might improve the DFT results close to the apex 
of the annular wedge. A more detailed investigation of tensor-weighted FMT with regard 
to the correlations in narrow planar slits and also in the annular wedges considered here is 
certainly of interest (although the algebra of App. [B] becomes much more extended in this 
case). 

Our observations should be compared with the phenomenological modelling of Ref. 26 1 
("annular slit approximation"). There, the effective 2d system was approximated by an 
idealized system of hard disks, layering around a hard cavity of radius tq centered around 
r|| = 0. The bulk density p2d,s of this hard disk system was determined via the self- 
consistency condition P2d(p2d,s) = — 27(ps) where p2d is the 2d pressure and 7 is the surface 
tension of a hard wall immersed in a bulk solvent of density ps- Employing scaled particle 



theory 



551], this yields pgds ~ 0-66 for p* = 0.764. This value is a bit smaller than the one 



inferred from the 2d density profiles of Fig. [61 The main difference between the idealized 



hard disk system of Ref. 



261 ] and the effective 2d system as reflected in Fig. [6] lies probably 



in the softness of the central cavity obstacle and also the softness of the effective particle 
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FIG. 6: (color online) Effective two-dimensional solvent density in the annular wedge for various 
solute-solvent size ratios. Lines correspond to DFT results (FMT-WBII) and symbols show results 
of our MC calculations. Left panel: solute-wall distance h = 0.95 a, right panel: h = a. The solvent 
bulk density is /o* = 0.764. The curves are plotted up to the point where the width of the annular 
wedge reaches the value of a. 



interactions due to the widening of the slit. Apart from that the overall picture of Ref. [26 1 
is confirmed well. In particular, this implies an important observable consequence: The 
cavity circumference (of approximate length 27rro) induces a line contribution 27rro 72d(p2d,s) 
to the insertion free energy of the solute near the wall. Since the insertion free energy is 
equivalent to the depletion potential up to an additive constant, the latter acquires a term 
oc ^^/bJo'^-Ii) , the corresponding term in the depletion force is oc \/R/ {a — h). (Note that 
for h ^ a the cavity radius should stay finite as reflected in Fig. [6] (right panel), thus the 
divergence of the depletion force there according to the geometrical argument of equating 
the cavity radius with tq is unphysical.) The interpretation of this line energy term within 
the more general framework of morphometric thermodynamics will be given below. 

We have seen that in radial direction, the packing of the solvent spheres close to the 
apex of the annular wedge is well understood by the quasi-2d picture developed above. For 
the particular point h = a, there is only moderate layering in r||-direction, thus the full 
density profile should be determined mainly by packing in ^-direction. Indeed, for h = a 
we observe a quite remarkable collapse of the density profiles at the wall when plotted as 
a function of f = r^^/y/Ra, see Fig. [71 indicating the relative unimportance of packing in 
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FIG. 7: (color online) Left panel: density at the wall for solute-wall distance h = a and the three 
solute-solvent size ratios 20, 30 and 100. Lines are FMT-WBII data, symbols are MC simulation 
results from Ref. [33(]. Right panel: The depletion force between solute and wall at h = a. The 



solvent bulk density is p* = 0.764. 

radial direction. According to Eq. ([2]), perfect scaling would imply that the depletion force 
at /i = 0" is proportional to R. The DFT data violate scaling for f < 0.2 (i.e. when the slit 
width is smaller than 0.02 a), as can be seen from Fig. [7] (left panel, inset). This leads to 
a slow increase of the scaled depletion force /u,/(2i?), see Fig. [7] (right panel). In contrast, 

indicate that /^/ (2i?) stays constant for large R, quite in 



the simulation results of Ref. 



agreement with an upper bound derived within the annular slit approximation of Ref. [26 1 
which reads /^/(2_R) = — 27r7(ps)(l— ?72d) where ?72d = (7r/4) p2d,sC"^ is the 2d "bulk" packing 
fraction. Note that the value of the Derjaguin approximation at h = a, /^/ (2i?) = — 27r7(ps) 
is completely off the simulation as well as the DFT values (see next subsection). 

2. Geometric interpretation of the depletion force 



Morphological (morphometric) thermodynamics 30|| provides a powerful interpretation 



of the depletion interaction, as has been shown very recently in the special case of the solute- 



solute interaction [29|]. Up to a constant, the depletion potential is nothing but the solvation 
free energy of the wall and the solute. In morphological thermodynamics the solvation free 
energy F^oi of a body is separated into geometric measures defined by its surface. As the 
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surface, we take the solvent-accessible surface, i.e. the body surface of the combined wall- 
solute object is defined by the dashed lines in Fig. [2] (a). These geometric measures are 
the enclosed volume V, the surface area A, the integrated mean and Gaussian curvatures C 
and X, respectively. To each measure there is an associated thermodynamic coefficient: the 
pressure p, the planar wall surface tension 7 and two bending rigidities k, and R such that 



Due to the separation of the solvation free energy into geometrical measures and geometry 
independent thermodynamic coefficients, it is possible to obtain the coefficients p, 7, k and 
R in simple geometries. The pressure p is a bulk quantity of the fluid. The surface tension 7 
accounts for the free energy cost of forming an inhomogeneous density distribution close to 
a planar wall. If the wall is curved the additional free energy cost is measured by k, and k. 
It is possible to obtain all four coefficients from a set of solvation free energies of a spherical 
particle with varying radius. For a hard-sphere solute in a hard-sphere solvent very accurate 
analytic expressions for the thermodynamic coefficients are known from FMT-WBII , see 
also App. |Al 

For the particular case considered here, the depletion potential W{h) follows by subtract- 
ing the sum of the solvation free energies of a single wall and a single solute sphere. Thus 
we find 



where AV^ and AA are the volume and surface area of the overlap of the exclusion (depletion) 
zones around wall and solute, respectively, and AC is the integrated mean curvature of that 
overlap volume. Note that the fourth characteristic, the integrated Gaussian curvature is 
the Euler characteristic which is An for wall and solute a.t h < a (one connected body) and 
Stt for wall and solute aX. h > a (two disconnected bodies). This explains the last term in 



i^soi = pV + -fA + kC + kX. 



(15) 



W{h <a) = -pAV - 'jAA - kAC - Atik 



(16) 
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FIG. 8: (color online) Depletion force between wall and solute for solvent bulk density p* = 0.764. 



Left panel: Solute-solvent size ratio 20. Simulation data are shown 
and differ by a slightly different choice of pw in Eq. ([2]) (see Refs. [31 



crosses (x) and circles (o) 



32] for details). Results from 



explicit FMT-WBII minimization are given by the full squares and results from lE-MV are shown 
by asterisks. The dash-dotted line shows the insertion route result using the WBII functional 
and the dashed line gives the Derjaguin approximation. The excess surface free energy of two 
parallel hard walls, needed for the Derjaguin approximation, has been calculated with the WBII 
functional. Right panel: depletion force for different solute-solvent size ratios. Symbols are results 
from explicit FMT-WBII minimization, lines in ascending order are morphometric results for size 
ratios 10, 20, 30, 50, 100 and oo (Derjaguin approximation). 



Eq. f[T6|) . The geometric measures AV, and AC are given by 
Ar = ^(a-/i)2(3i?-(a-/i)) 

""-"^ TxR{a-hf + 0{R'') , 
AA = A-KR{a -h)- Tx{a - hf 

"^^^ 47ri?(fT -h) + 0{R'') , (17) 
AC = 2n{a - h) + n^/2R{a - h) - (a - hf + arcsin - ^ 

"""" 7rV2i?(a-/i) + 0(i?°) . 



The depletion force follows as fw{h) = —dW/dh. In the limit i? — oo it is given by: 



2nR ' ' 2 V 2R{a - h) ^ ' 
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The contribution from the integrated mean curvature is subleading in but it is quantita- 
tively important since it induces a singularity as h ^ a which is proportional to l/(cr — /i)^/^. 
This contribution precisely corresponds to the line tension term associated with the circular 
apex (of length 27rro) of the annular wedge identified in the quasi-2d analysis of Ref. 26 1 
and introduced in the last subsection. The line tension in the quasi-2d picture was shown 
to be 72d(p2d,s) and in the geometric analysis it is —7tk/2. There is good agreement between 



both expressions 



29|. Strictly, the singularity associated with this line term is unphysical 



of course. Indeed, it can be shown that the morphometric analysis becomes ambiguous very 
close to h = a. One would require that the morphometric result for the force is unchanged 
if another, slightly displaced surface around the solute sphere is chosen for the geometric 



analysis. In Ref. 



291], such an analysis was carried out for the sphere-sphere geometry and 



the independence of the force on the chosen surface was demonstrated. However, if one 
chooses e.g. the molecular surface (the surface of the set of points which is never covered 
by a small solvent sphere), then this surface becomes self-overlapping close to h = a and 
one could not expect to ascribe physical significance to surface tensions and mean curvature 
coefficients of such overlapping surfaces. For the sphere-wall geometry a similar analysis 
holds. Note that the singularity at h = a corresponds to a vanishing circumference of the 
apex (ro = 0), however, the 2d analysis of the density distributions at this point yielded a 
small but finite value for the "effective" apex radius ro (see Fig. [6] (right panel)). 

In Fig. [8] (left panel) we show FMT-WBII and lE-MV results for the depletion force 
fw between a big sphere and a wall (size ratio 20, solvent bulk density p* = 0.764) and 
compare them to MC results from Refs. 3l|, l32|] as well as to the morphometric analysis. 
The force according to morphometry is plotted only until self-intersection of the molecular 
surface starts. The FMT data are described very well by the morphometric results; the sharp 
drop in the depletion force close to h = a appears to mimick the behavior of the singular 
term oc dAC/dh. The MC data suffer from relatively large error bars (except for the point 
h = a but are overall consistent with both the FMT data and the morphometric analysis. 
The scatter in the lE-MV data is due to the comparatively low resolution in 2;-direction 
as compared to the DFT data {Az = 0.05 a vs. Az = 0.002 cr). As before, all considered 
IE closures give similar results, but these correspond to an overall more attractive force 
compared to DFT and MC. Two other approximations for the depletion force are shown 
in Fig. [8] (left panel): the DFT insertion route and the Derjaguin approximation. Both 
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approximations fail to describe the depletion force near h = a. In the insertion route to 
DFT (see Ref. for its derivation), one circumvents the explicit calculation of the density 
profile around the fixed wall and solute by exploiting the relation: 



/3W{x.) = lim 



-f^TiPs). (19) 

Ps{r)=Pw{z) 



Here, pl^{ps) is the excess chemical potential for inserting one solute into the solvent with 
density ps- The functional derivative of the mixture functional has to be evaluated in the 
dilute limit for the solutes, pb ^ 0, and with the solvent density profile given by the profile 
at a hard wall {py^{z)). This method relies on an accurate representation of mixture effects 
for large size ratios in the free energy functional. Therefore, the observed deviations in the 
depletion force presumably follow from the insufficient representation of higher-order direct 
correlation functions of strongly asymmetric mixtures in the present forms of FMT 57|]. The 
Derjaguin approximation, on the other hand, corresponds to the leading term for i? — >■ 00 



in Eq. ffTSj) inside the depletion region {h < a), but is easily extended to h > a j25|, l26| : 

jDerjaguin J ^(/^ _ o") - 27 {h < a) 



27tR 



(20) 

7siit(/i) - 27 {h> a) 



with 7siit(^) denoting the excess surface energy of two parallel hard walls at distance h 
forming a slit. According to Eq. ([21), the Derjaguin force scales with R ("colloidal limit") 
which is an excellent approximation outside the depletion region. (Incidentally, outside the 
depletion region all methods and the Derjaguin approximation agree with each other.) Inside 
the depletion region, the neglected mean curvature (or line tension) term is very significant. 

Explicit FMT-WBII results for the depletion force for size ratios up to 100 are shown 
in Fig. [8] (right panel) and compared to the morphometric analysis (symbols vs. full lines). 
Whereas for size ratio 10 some discrepancies are still visible, they have more or less vanished 
for size ratio 100 (except for the point h = a, see the previous discussion). Thus the 
morphometric form for the depletion force (Eq. (fTSl) ) can be regarded as the appropriate 
"colloidal limit" . 



B. Solute solute interaction 



We now turn to the case of the depletion force between two big solute particles which are 
immersed in the hard solvent. Considering the evidence from the sphere- wall case, we would 
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FIG. 9: (color online) Left panel: Depletion force between two solutes. Circles correspond to 
the centered difference of the FMT-WBII data for the depletion potential (right panel), while 
squares give the FMT-WBII depletion force according to the surface integral in Eq. ([1]). Right 
panel: Depletion potential between two solutes. The curve for the morphometric potential has 
been shifted by -1. The solvent density is p* = 0.8 and the solute-solvent size ratio is 5. 

expect that for large size ratios, the morphometric picture reliably describes the depletion 
force. Since the geometry of the overlap volume between two spheres is different from that 
between a sphere and a wall, we find also a slightly different form for the morphometric 
force in the limit R ^ oo: 



f{h<a) R^oo X „ TV I 

ii-^ « p(A^,)_2o^_,_^__. (21) 

Thereby, one sees that the scaling relation Z^^, = 2/, valid within the Derjaguin approxi- 
mation, is violated through the appearance of the mean curvature (line tension) term. For 
moderate solvent densities up to p* = 0.7 the morphometric form ( 12T|) gives a good descrip- 



tion of the depletion force for size ratios a > 10. This has been reported in Ref. 29|]. An 
analysis of the full solvent density profiles between the two solutes reveals the same features 
as described in Sec. IIII A[ 

It is interesting to investigate the regime where the morphometric description is expected 
to fail. This will happen when there are strong correlations in the whole annular wedge, 
i.e. where packing in both the r\\- and the 2;-direction becomes important. Typically, 
intermediate size ratios and higher solvent densities will induce these strong correlations. 
Therefore, we have investigated the depletion force for size ratios 5 and 10 and solvent 
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densities p* = 0.8 and p* = 0.9. In Figs. [9] and [TO] we show results for the depletion force 
and potential at p* = 0.8 for size ratios 5 and 10 respectively. For size ratio 5, the explicit 
DFT calculations yield a depletion force which does not correspond to the morphometric 
result at all. It is closer to the insertion route (see (Eq. [I9])), though systematically lower. 
This yields a depletion potential (Fig. [9] (right panel)) which is about 1 k-sT or 25% more 
attractive at contact than that of the insertion route. For size ratio 10, the explicit DFT data 
for the depletion force reflect the morphometric form a.s h a, i.e. they give again evidence 
for importance of the mean curvature (line tension) term in Eq. (I2T!) . However, away from 
that regime, the depletion force deviates significantly from both the morphometric form 
and the insertion route result such that the depletion potential at contact (Fig. (right 
panel)) is about 3 k-sT or 50% more attractive at contact than that of the insertion route 
(for FMT-WBII). 

Note that at such a high solvent density (p* = 0.8) there are already significant deviations 
between FMT-RF and FMT-WBII (Fig. [10] (right panel), circles and squares). Since FMT- 
WBII has been designed to improve thermodynamic and morphometric consistency at higher 
densities, the corresponding results can be assumed to be more trustworthy. As an additional 
check of the numerics, we calculated the depletion force in two ways: (a) via the surface 
integral in Eq. ([1]) and (b) via the centered difference of the results for the depletion potential 
(see Figs. [9] and [10] (left panel), squares and circles). The latter is simply obtained as 

W{h) = fi[p(r); /^]|^(r)=p.,(,„.;.) - mry,h - oo]|^(r)=p.,(.,„.;.^oo) , (22) 

where the grand potential Q and the equilibrium density peq(r||, z; h) around the two solutes 
at distance h are determined by the basic DFT equations ([3]) and ([5]), respectively. The 
agreement between routes (a) and (b) is very good, save for some "fluctuations" in route (b) 
close to the point h = a where the density oscillations in the wedge are most pronounced. 
The equality of both routes checks a sum rule check similar to the hard wall sum rule for the 

problem of one wall immersed in solvent. Weighted-density DFT's usually fulfill the latter 

58|. 

The differences between the insertion route and the explicit FMT data becomes more 
dramatic for higher densities. In Fig. [TT] we show results for the depletion potential for a 
solute-solvent size ratio 5 at solvent densities 0.85 and 0.9 (left panel) and for a size ratio 
10 at density 0.85 (right panel). For size ratio 5 the repulsive barrier has almost vanished 
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FIG. 10: (color online) Left panel: Depletion force between two solutes. Circles correspond to the 
centered difference of the FMT-WBII data for the depletion potential (right panel), while squares 
give the FMT-WBII depletion force according to the surface integral in Eq. ([1]). Right panel: 
Depletion potential between two solutes. Squares correspond to FMT-WBII results, and circles to 
FMT-RF results. The curve for the morphometric potential has been shifted by -2. The solvent 
density is p* = 0.8 and the solute-solvent size ratio is 10. 

for p* = 0.9, and the potential well at contact is almost twice as deep as compared with the 
insertion route. For size ratio 10, the barrier remains (albeit at a different location) but the 
well depth is equally enhanced by a factor 2 compared with the insertion route as for size 
ratio 5. 

The significant deviations of the depletion potential (according to the explicit FMT- 
WBII results) from the previously known approximations (Derjaguin approximation and 
insertion route) sheds some doubts on the quantitative accuracy of simulations of binary 
hard sphere mixtures (with asymmetries between 5 and 10), employing an effective one- 
component approach for the solutes interacting with their depletion potential. In Ref. |l^ 
the phase diagram of the mixture was scanned in that way using a virial expansion to 
third order in the Derjaguin approximation for the depletion potential (see the dot-dashed 
curves in Figs. [9] and [10] (right panel)). Especially for size ratio 5, the explicit DFT results 
predict a depletion potential with much smaller barrier and a deeper well at contact. This 
might affect the phase diagram in the corner where the packing fraction of the solutes 



is low and the solvent packing fraction is high. In Ref. 59|] gellation in hard sphere-like. 



asymmetric mixtures is investigated with an integral equation approach (reference functional 
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FIG. 11: (color online) Depletion potential between two solutes, comparison between explicit FMT- 
WBII data and the insertion route. Left panel: Solute-solvent size ratio 5 and solvent densities 
p* = 0.85 and 0.9. Right panel: Size ratio 10 and p* = 0.85. 



approximation 



ll| ) to the depletion potential, akin to the insertion route. We expect 



quantitative deviations at higher solvent densities compared to explicit DFT calculations 
which again affect the corner of the phase diagram with low solute and high solvent packing 
fraction. 

Finally we want to mention that the explicit DFT methods introduced here, together with 



the morphometric 



analytic forms 
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'orm of the depletion force in Eq. ( |2T]) . can be used to improve available 



6l| for the contact values of the pair correlation functions in hard-sphere 



mixtures (note that the solute-solute contact value is given by gbbi^2Rb) = exp[— /3Vr(0)]). 



IV. SUMMARY AND CONCLUSIONS 



We have presented a detailed analysis of the relationship between depletion forces (solute- 
wall and solute-solute) in a hard solvent and the associated density distributions around 
the fixed wall/solute particles, using mainly density functional theory of fundamental mea- 
sure type complemented with Monte-Carlo simulations and integral equation techniques. 
For large solute-solvent size ratios, the depletion force is strongly linked to the quasi two- 
dimensional confinement between the solutes (or the solute and the wall). We have shown 
that the properties of this quasi two-dimensional confined system are reflected in the de- 
pletion force by a line-tension term which in turn can be obtained quite generally through 
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morphological thermodynamics. Thus we can formulate an appropriate "colloidal" limit for 
the depletion force in hard-sphere mixtures: Eq. f|T8l) for the solute-wall case and Eq. fl2T|) 
for the solute-solute case. This new formulation improves significantly over the Derjaguin 
approximation which is a frequently employed tool to estimate colloidal interactions in many 
circumstances. 

Although our analysis has been restricted to hard spheres only, the formulation of mor- 
phological thermodynamics is very general such that it can be expected that the colloidal 
limit for the depletion force holds also in mixtures with more general interparticle interac- 
tions. However, more explicit tests in this direction are necessary. 

For intermediate size ratios 5 and 10 and higher solvent densities p* > 0.8 we found strong 
deviations in the solute-solute depletion force between the explicit density functional results 
on the one hand and the Derjaguin approximation or morphometry on the other hand. This 
is a result of the strong solvent correlations in the annular wedge between the solutes. The 
depletion well at solute contact is significantly more attractive (by almost a factor of two) 
than previously estimated. This might have consequences for the phase diagram of binary 
hard spheres at low solute and high solvent packing fractions. With recently developed 



techniques 
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63| . this question appears to be tractable in direct simulations. 
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APPENDIX A: FUNDAMENTAL MEASURE FUNCTIONALS FOR HARD 
SPHERES AND ASSOCIATED THERMODYNAMIC COEFFICIENTS 

We consider fundamental measure functionals involving no tensor-weighted densities 
which are defined by the following functional for the excess free energy of a hard sphere 
mixture of solvent and solute particles, described by the density distribution p(r) = 
{Ps{r),Pb{r)}: 

(ir- = I cir$({n[p(r)]}), (Al) 

*({n[p(r)]}) = -no ln(l - ris) + (^1(713) + <^2{n3) ^. ,^ ^ ■ 

1 — 77.3 /47r(l — riz) 
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Here, $, is a free energy density which is a function of a set of weighted densities {n(r)} = 
{no, ni, 77-2, ^3, 111, with four scalar and two vector densities. These are related to the 
density profile p(r) by 



J2 I dv'p.iv') wtir -r') = J2p,* < , (A2) 



and the hereby introduced weight functions, {w,j(r)} = {Wi,wl,wf,wf,w},wf}, depend on 
the hard sphere radii Ri = {Rg, Rb} of the solvent and solute particles as follows: 

wf = eiR,-\v\), w^ = 5{R,-\v\), w] = ^, w'- 



w,^ = ^MR. - |r|) , w,H . (A3) 

|r| AnRi 

The excess free energy functional in Eq. flAll) is completed upon specification of the functions 
(pii^ri^) and (p2{ns). With the choice 

(^1 = 1, (^2 = 1 (A4) 

we obtain the original Rosenfeld functional j^. Upon setting 

^1 = 1 (A5) 
-2n3 + 3n2- 2(1- 723)2 ln(l-n3) 

(P2 = I 

we obtain the White Bear functional 
Starling equation of state. Finally, with 



3ni 



5|, consistent with the quasi-exact Carnahan- 



2n3-nl + 2(l-n3)ln(l-n3) 

ipi = l-\ (A6) 

3n3 

2r23 - 3nl + 2nl + 2(1 - ^3)^ ln(l - ^3) 

"'^^ 

the recently introduced White Bear II functional is recovered. 

Next we briefly recapitulate the determination of the thermodynamic coefficients p (pres- 
sure), 7 (hard wall surface tension), k and R (bending rigidity of the integrated mean and 
Gaussian curvature, respectively) |6| . Consider the free energy of insertion Fi, for a big solute 
particle of radius Rb and associated radius R = Rb + Rs of the exclusion sphere around it. 
According to Eq. (fTSll . it is given by: 

47r 

Fb = p—R^ + -fATrR^ + KATrR + RATr. (AT) 
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On the other hand, Fi, is equivalent to the excess chemical potential of the solutes in a 
mixture with solvent particles at infinite dilution {pb 0). This excess chemical potential 
is obtained as the density derivative d/dpb of the mixture free energy (Eq. (lAip ). therefore 
we obtain: 

By equating the two expressions for Fi, in eqs. (]A7p and (lASp . one finds the explicit expres- 
sions for the thermodynamic coefficients as linear combinations of the partial derivatives 
d^/dui. For the Rosenfeld functional, these coefficients are given by: 

1+Vs+ Vs 



PP = Ps 



;i - Vs 



^3 



3 r]s{l + r]s) 



(3k = pscr 



4"^ {1-Vsr 

2 n 

—p«cr 

4^' (1 - r^s 

2 + 7r)s- llry^ ln(l - r/. 



48(1 - r/,)3 24r/, 

For the White Bear II functional, the thermodynamic coefficients have been given already 
6| and are reproduced here for completeness: 

I + r]s + ril - 4 



in Ref. 



Pp 



p-f = -psCr 



1 - Vs)' ' 

l + 2r], + Srjl - hvil ln(l - r/. 



6(1 - risY QVs 



Pn = p^a- ( ' '-^^-/^Qf + ) , (AlO) 



/3k = psO 



6(1 - risY 'iri. 
-4 + ll?7s - 13?7^ + 4?7f ln(l-?7, 



24(l-r/,)3 6r/, 

The coefficients of the White Bear II functional are remarkably consistent with respect to 
an explicit minimization of the functional around solutes of varying radius and fitting the 
corresponding insertion free energy to Eq. (1A7P 6|. 
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APPENDIX B: MINIMIZATION OF FUNDAMENTAL MEASURE FUNCTION- 
ALS IN CYLINDRICAL COORDINATES 



In the following we consider a one-component fluid of solvent hard spheres only (ps(r) = 
p(r)). The equilibrium density profile Pcq(r) of the hard sphere fluid with chemical potential 
fi = (3~^ ln(ps A'^) + (corresponding to the bulk density ps) in the presence of an arbitrary 
external potential V^(r) is found by minimizing the grand potential 

m = -^'"[P] + -^'^^[P] - / dr{p - V^^Xv)) , (Bl) 



which leads to 



^-1 1^ PcqW ^ + p<=^ - K^^*(r) . (B2) 

Ps 

The functional /i[p(r)] is given by 

- ^ (B3) 

For the physical problems of this work, sphere-sphere and wall-sphere geometry, the external 
potential possesses rotational symmetry around the z-axis. Therefore we work in 
cylindrical coordinates r = (r|| cos 0, r|| sin 0, z), in which the external potential and the 
density profile depend only on r\\ and z, V = V{r\\,z) and pcq = Pcq{r\\,z). Eqs. flB2p and 
(]B4p can be solved with a standard Picard iteration procedure or a speed-enhanced scheme 
(see below) . The technical difficulty lies in the fast and efficient numerical evaluation of the 
weighted densities Ua = p * w"' and the convolutions d^/dua * w"" appearing in Eq. (1B4I) . 

Convolution integrals are calculated most conveniently in Fourier space where they re- 
duce to simple products. Generically, the (three-dimensional, 3d) Fourier transforms of the 
weighted densities involve a one-dimensional (Id) Fourier transform in the 2;-coordinate and 
a Hankel transform of zeroth or first order in the r||-coordinate (see Sec. IB ip . Both the Id 
Fourier transform and the Hankel transform can be calculated using fast Fourier techniques 
(see Sec.|B2]). 
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1. Weighted densities 



The convolution integral, defining the weighted densities in Eq. (]A2p . is calculated by: 

j^--^exp{-iq ■ r) p{q\\,qz) w(q) . (B5) 

The 3d Fourier transform of the density profile p(r|[, z) appearing in the above equation is 
given by: 

p{%,Qz) = j cirexp(iq- r)p(r||,2;) (B6) 

/OO noo 

dzexp{iqzz) / 27rr||dri| Jo(gi|r||) p(r||, 2;) , (B7) 

-00 Jo 

= FTHTop(r||,z) . (B8) 

Here, we have introduced shorthand notations FT for the Fourier transform in the z- 
coordinate and HTj for the Hankel transform in the r||-coordinate involving the kernel 
Ji{q\\r\\). A vector in real space is given by r = r||e|| + ze^, whereas a vector in Fourier space 
is given by q = gne'i + g^e^, with e|| ■ ej| = cos0g. The 3d Fourier transforms w of the set 
of weight functions w reduce to sine transforms due to radial symmetry and are explicitly 
given by 



-^(^) = ^ - cos(g/2) j , (B9) 

w\q) = ^ sm{qR) , (BIO) 
w2(q) = -iqw^q) (Bll) 

(B12) 

(The remaining weighted densities differ only by a multiplicative factor.) It is convenient to 
introduce the following parallel and perpendicular components of the Fourier transformed 
vector weights {k = 1,2): 

= ^e|| ■ w*^(q) , (B13) 
= le, ■ w'=(q) . (B14) 

Using these definitions, the scalar weighted densities are given by 

n,(r||,^) = FT"iHTo"^[p(g||,g.)w}'=(y/^^T^)] . (A; = . . . 3) . (B15) 
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The vector weighted densities, on the other hand, are given by two components {k = 1, 2): 

nfc(r) = nk,\\{r\\, z)e\\ + nk,^{r\\, z)e^ , 

nk,\\{r\\,z) = -FT-iHTi-i[p(g||,g,)*[f(g||,g,)] , (B16) 

nk,zi^\\,z) = FT "^HTo . 

We demonstrate this resuh for the example of the weighted density n2. Upon choosing 
r = (r||, 0, z) we find 



n2(r) 



(iq 

(27r)3 

^ dq, 



exp(-zq • r) p(gii, z) (-zq)w; (g) 



(BIT) 



(2vr) 



- exp{-iq^z) I q\\dq\\p{q\\,z) j rf^g exp(-zg||r|| cos </)g) 



2-K 



iw"^^ cos< 








sin <; 


K 






-iwl 


J 



^ -Ji{q\\r\\)w'^^ 




(B18) 



\^-iJo(g||r||)w2y 

which is equivalent to Eq. flB16p for k = 2. Here we made use of the integrals 

/ (i0exp(— COS 0) = 27rJo(x) , 
Jo 

p2TT 

/ (i</)exp(— ix cos</)) sin</) = 0, 
Jo 

r-2-K 

/ (i0exp(— ZXCOS0) cos</) = ^-nU'^ix) = —2mJi{x) . 
Jo 

Next we consider the evaluation of the convolution type integrals appearing in fi[p\ (see 
Eq. dEH)): 

/^M = 5Z / dr'p^{r')w"{r' - r) , (B19) 

where pa = d^/dua and a runs over scalar and vector indices. For scalar indices, w^{r) = 
w^{—r): we recover the standard convolution integral and thus: 

(A; = . . . 3) . (B20) 



dT'pk{r')w\r' - r) = FT "^HTo [pk{q\\,q^) [^q^ + ql 

In the case of vector indices, we observe that the free energy density $ only depends on 
vector densities through rii ■ n2 and n2 ■ n2. Thus we find {k = 1, 2): 

Pk{T^) = Pk,\\{r\\,z)eii+pk,zirii,z)e^ (B21) 

Pfc(q) = Pk,\\{q\\,qz)e'n+pk,z{q\\,qz)ez (B22) 
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with 

Pk,\\{q\\,qz) = iFT HTi pk,\\{r\\,z) , Pk,z{q\\,qz) = FT HTq pk,z{r\\,z) . (B23) 
Therefore the vector part summands of fi[p] are evaluated by 

J dr'pk{r')w\r' - r) = FT -^HTq [-pkM^ + zpk,zW^] (B24) 

2. Fast Hankel transforms 

Hankel transforms can be calculated by employing fast Fourier transforms on a logarith- 
mic grid. Consider the Hankel transform 

/"OO 

HT^/(r||) = 27T / riidrii J^(g||r||)/(r||) . (B25) 
Jo 

We define new variables x = ln(r||/ro) and y = ln(g||/go) with arbitrary constants Tq and go- 
In terms of these variables 



oo 



HT^ /(rii) = 27rro' / dx J,{x + y) f{x) , (B26) 

J — oo 

where J^(x) = J^{qoroe^) and /(x) = e^^fi^e^). Thus, the Hankel transform takes the 
appearance of a cross-correlation integral and can be solved via Fourier transforms: 

HT^ /(rii) = 27Trl FT [fT J^{x) FT*/(x)] . (B27) 

The Fourier transform of can be done analytically while the remaining ones are computed 
numerically using Fast Fourier techniques. 



In the actual implementation we followed Ref. [6J] which details the proof of orthonor- 
mality for discrete functions, defined on a finite interval in logarithmic space (and continued 
periodically). The following remarks apply: 

• The fast Hankel transform cannot be applied directly to the density profile, f{r\\) = 
p{r\\, z), since it does not go to zero for r\\ oo. It is therefore advantageous to split 
the external potential 

^-t(r) = V^(^) + K(r||,^) (B28) 
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into a (possibly) ^-dependent background part V{z) and the remainder. The corre- 
sponding background profile pv fulfills 

In = -f3^[py(z)] + _ f3V{z) . (B29) 

ps 

The full equilibrium profile can then be written as pc(i{r\\, z) = pv{z){h{r\\, z) + 1) with 
h{r\\,z) for r\\ oo. The determining equation for h reads 

Hh{r\\,z) + 1) = -I3p[pv{z)h{r\\,z) + pv{z)] + I3p[pv{z)] - [3V,{z) . (B30) 

Thus one needs to perform Hankel transforms only on functions which properly go to 
zero for r\\ —>■ oo. In the wall-sphere case, V{z) is naturally given by the wall potential, 
Vs becomes the solute-solvent pair potential Ubs and h = hbs is the solute-solvent pair 
correlation function. In that form, Eq. (]B30p resembles the general closure for integral 
equations (Eq. ffTOj) ). In the sphere-sphere case, V{z) = with pv = Ps- 

The assumed periodicity of f{x) leads to restrictions on the product go'^o (low-ringing 



condition in Ref. 6J]). This condition cannot be fulfilled on one grid for both HTq (■) 
and HTi (■). However, this does not lead to any noticeable instabilities in the repeated 
application of the fast Hankel transform. 

In order to avoid aliasing and the amplification of numerical "noise" in the low-x and 
high-x tails of f{x) in the repeated application of the fast Hankel transform we worked 
with cutoffs rmin/ro = gmin/go = 0.01 and r^ax/^o = ^max/go = O(IOO). The fast 
Hankel transform itself was calculated on an extended grid x G [—A^|| Ax/2, A^|| Ax/2] 
with either A^n = 2048, Ax = 0.01 or A^n = 4096, Ax = 0.005. Outside the "physical" 
domain defined by Xphys G [ln(rmin/ro), ln(rmax/'ro)] the function /(x) was put to zero 
(a similar prescription applies to 2/phys)- 



3. Speedup of iterations 

The density profile /o(t'||, z), the weighted densities n(r||, z) and the derivatives p(?"||, z) = 
d^/dn{r\\, z) have been discretized on a two-dimensional grid spanning the plane {r\\,z). 
Spacing in ^-direction was equidistant with grid width Az = 0.002 . . . O.OOScx with up 
to Nz = 15,000 points. Spacing in ry-direction was logarithmic (see above) with N\\ ^ 
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1, 000 ... 2, 000 points. Memory requirement went up to 16 GB for the largest grids. Com- 
putations have been performed on nodes with 16 GB RAM and two Intel QuadCore pro- 
cessors, with OpenMP parallelization of the arrays of Fourier and Hankel transforms. One 
evaluation of yu[p] (one iteration) took up to two minutes. 

The equilibrium density profile fulfilling Eq. (lB2p or Eqs. (1B29P and (1B30P can be de- 
termined by Picard iterations where, as a minimum requirement for convergence, careful 
mixing of the current and previous iteration is necessary. However, for packing fractions 
7] > 0.3 and larger colloid-solvent size ratios a > 5 one easily needs several hundreds to 
thousands of iterations until convergence. In view of the iteration times up to two minutes, 
this is inacceptable. Therefore we ernployed the modified DIIS (direct inversion in iterative 
subspace) scheme developed in Ref. 65| which essentially constructs the next iterative step 
out of a certain number of previous steps. In the DIIS scheme, the mixing coefficients of 
the previous steps, determining the solution of the next step, are obtained by a minimiza- 
tion condition on the residual. The modification to DIIS consists in the admixture of the 
extrapolated DIIS residual to the solution of the next iterative step, in order to enlarge 
the dimensionality of the iterative subspace and thus to reach the true solution much more 
quickly. Our practical experience with modified DIIS is very similar to the observations in 



Ref. 



65|, this includes the necessity to combine DIIS with Picard steps carefully, in case the 



DIIS steps show divergent behaviour. In summary, modified DIIS is a robust method for 
our problem, and the total number of iterations is reduced to approximately 100 and less for 
most parameter choices. The only noticeable exception occured for the case of the nearly 
singular annular wedge, i.e. when the distance h between wall and colloidal sphere (or the 
two colloidal spheres) is ~ 1 a and the solvent packing fraction 77^ > 0.35. Here, several 
hundred of intermediate Picard steps between two DIIS steps where occasionally necessary 
to prepare the ground for the next DIIS step. 
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APPENDIX C: INTEGRAL EQUATIONS FOR THE SOLUTE WALL CASE: 
NUMERICAL SOLUTION 



In cylindrical coordinates the inhomogeneous OZ equations ([7j)-([H]) read 

hss{z,zo,s) - Css{z,zo,s) = j ds'dz'pv{z')hss{z,z',\s-s'\)cssiz',zo,s') , (CI) 

hhs{z,zo,s) - Chs{z,Zo,s) = J ds'dz'pviz')hbs{z,z',\s- s'\)css{z',Zo,s') , (C2) 

where the total and direct correlation functions hij{z, zq, s) and Cij{z, zq, s) of two particles 
depend on the distances z and zq from the wall and on the projection of their position vectors 
on the direction along the wall, s = r|| — r|| q- The integral over s' is a 2D convolution, which 
is most efficiently done by means of the Fast Hankel transform of the zeroth-order (see 
Appendix IB 21) . In Fourier space, the OZ equations become 

hij{z,Zo,q\\) - Cij{z,Zo,q\\) = ^ / dz' pv,k{z')hik{z, z' ,q\\)ckj{z' , Zo,q\\) , (C3) 

k=b,s 

where h[c]{z, zo,q\\) = HTq h[c]{z, zo,r\\). The remaining ^-integral is evaluated with a 
simple trapezoidal rule on a uniformly discretized grid of Nz points, yielding the following 
matrix equation: 

H,,(g||) - C.,(g||) = J2 H.fc(g||)RC,,(g||) , (C4) 

k=b,s 

where H and C are A'^^ x A'^^ matrices generated by the corresponding correlation functions 
at each g|| point and R is a diagonal matrix corresponding to py{z)6{z — z'). 

The Lovett-Mou-Buff-Wertheim (LMBW) equation (Q for the background density pro- 
file in the presence of a hard wall takes the following form: 

^^ = py(^) [pv{0) I rfr||c..(^,0,r||) + J dz'dv^^cUz, z' ,r\\)^P^^ . (C5) 

Here the 2D integration over r\\ can be considered as the DC component q\\ = of the 
Hankel transform of the direct correlation function Css{z , z' , q\\) and the remaining integral 
in the second term is taken again with the trapezoidal rule. Note, that the contact density 
at the plannar wall is known exactly from the wall theorem pv'(O) = Pp, where p is the bulk 



pressure 66|]. Thus, a substitution of the contact density by the expression for the bulk 



pressure provides the modified LMBW equation for the density profile 67|. Following the 
ansatz of Plischke and Henderson 68|], we used the Carnahan-Starling fit for pressure, which 
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is the first equation in (lAlOp . This in general should improve the accuracy of the density 
profiles near a plannar wall, if an approximation to the solvent direct correlation function 
Css{z, z' ,r\\) is applied. 

The closure equation flTU]) is the central one in the elaboration of a successful theory, 
since the OZ and LMBW equations described above are formally exact and need to be 
complemented with a third relation between between correlation total and direct functions. 
Unfortunately, in its exact form this third relation, expressed via additional bridge function 
6jj(r, fq) is defined only as an infinite sum of highly connected bridge diagrams and so cannot 
be fully utilized. In practice, one needs to resort to different approximations, which are suited 
for particular applications depending on the system state. We considered the Percus-Yevick 
(PY), Rogers- Young (RY) and modified Verlet (MV) closures defined in Eqs. ffTTl) - f[T5]) . For 
RY we use a scaled one-parameter form for ^jj, = ^/ {Ri + Rj) with ^ = 0.160 to fulfill the 
single-component thermodynamic consistency requirement. For ^ = one recovers the PY 
closure (Eq. lfTTl) ). while for ^ — oo one obtains the HNC closure. Likewise the HNC closure 
is recovered if — > oo in the MV closure (Eq. ( fT3l) ). MV is consistent with PY up to the 
fourth virial coefficient, if the aij are density independent. Here we follow the suggestion 



of Henderson et. al. 52| and use their definition for the state-dependent parameters aij. 
In the infinite dilution limit considered in the present paper these parameters reduce to the 
one-component form 

a = exp(-2''= ) + 0.8 - 0.45r/, . 

120?7s 

The contact values from the MV closure proved to be in good agreement with the val- 
ues calculated from Carnahan-Starling equation of state, which constitutes a considerable 
improvement over the PY and HNC approximations. 

In the application of these closures to the inhomogeneous system the following precautions 
must be observed to obtain robust results. 

• Technically we use the same iteration scheme as in Sec. IB 31 for the numerical solution 
of Eqs. flC4l) - flC5l) and ffTOl) (though due to the memory restrictions the z-grid spacing 
was only Az = 0.05a with up to A^^ = 400 points). The convergence drastically 
depends on the initial conditions, i.e. the zeroth iteration. While solving the PY 
closure it was enough to initialize both solvent-solvent and solute-solvent correlation 
functions with zeros, the RY and MV solute-solvent correlations have to be initialized 
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with a good guess to the final solution (which was actually the corresponding PY 
result). The LMBW equation imposes even stronger requirement on the initial profile, 
which needs to be either close to the true profile or evolve very slowly from the flat 
profile together with the corresponding correlation functions. 

• It turned out that the order of iterations is crucial for the overall convergency: for a 
given approximation to the density profile one needs to solve the closure and OZ equa- 
tions prior to the next iteration on the LMBW equation. It is therefore very costly to 
get a fully self-consistent solution of all three equations and normally the convergence 
goal for the background density pv{z) is much lower than for the correlation functions. 

• To avoid the necessity of requiring the bulk equilibrium behavior in correlation func- 
tions for distances far enough from the hard wall, we introduced the second wall at 
L ~ 20(T apart from the first one. The wide slit geometry allowed us to keep the z- 
resolution at the maximal possible, yet feasible level. Finer resolution in z in general 
improves numerical stability of the iteration procedure. Alternatively, one can use a 
very recent method of the expansion into the orthogonal set of functions proposed by 
Lado in Ref. M . 
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